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Abstract 

We study the quantum many-body dynamics and the entropy production triggered by an inter¬ 
action quench in a system of N = 10 interacting identical bosons in an external one-dimensional 
harmonic trap. The multiconfigurational time-dependent Hartree method for bosons (MCTDHB) 
is used for solving the time-dependent Schrodinger equation at a high level of accuracy. We con¬ 
sider many-body entropy measures such as the Shannon information entropy, number of principal 
components, and occupation entropy that are computed from the time-dependent many-body ba¬ 
sis set used in MCTDHB. These measures quantify relevant physical features such as irregular or 
chaotic dynamics, statistical relaxation and thermalization. We monitor the entropy measures as 
a function of time and assess how they depend on the interaction strength. For larger interaction 
strength, the many-body information entropy approaches the value predicted for the Gaussian 
orthogonal ensemble of random matrices and implies statistical relaxation. The basis states of 
MCTDHB are explicitly time-dependent and optimized by the variational principle in a way that 
minimizes the number of significantly contributing ones. It is therefore a non-trivial fact that sta¬ 
tistical relaxation prevails in MCTDHB computations. Moreover, we demonstrate a fundamental 
connection between the production of entropy, the build-up of correlations and loss of coherence 
in the system. Since the coherence and correlations are experimentally accessible, their present 
connection to many-body entropies can be scrutinized to detect statistical relaxation. Our re¬ 
sults are the first ones obtained for thermalization of finite quantum systems using an optimized 
time-dependent and genuinely many-body basis set. 

PACS numbers: 05.45.Mt, 02.30.1k, 05.30.-d, 05.70.Ln 
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I. INTRODUCTION 


The onset of thermalization in an isolated quantum system with a finite number of in¬ 


teracting particles is an important issue in theoretical physics in recent days 


15j. Ex¬ 


perimental progress with various interacting quantum systems has further corroborated the 


interest 


16 


171 ]. These experiments reveal that the question of how entropy is produced 


m 


a quantum system is one of the basic outstanding problems of many-body physics CHID- 
A necessary condition for thermalization is statistical relaxation in various observables of 
a system to some kind of equilibrium 

1 Il2i| 


interacting particles, statistical relaxation is related to chaos in the energy spectra 


20 


and here the inter-particle interaction plays an important role. Chaos in quantum systems 
is in turn defined by the statistics of eigenstates. For time reversal and rotationally invari¬ 
ant systems these follow the predictions of the Gaussian orthogonal ensemble of random 
matrices (GOE), see Refs. 2, 3, 18, 21]. 

The eigenstate thermalization hypothesis (ETH, put forward in Ref. j^j]) was so far only 
tested using models that employ a fixed basis set related to noninteracting particles. For 
example, calculations have been done for one-dimensional spin- 1 systems with nearest- 


neighbor anc 

n n 

bosons 9, 18 


next-nearest-neighbor coupling as well as for gapped systems of hard-core 


2l| . Statistical relaxation was found, for large enough interaction strength, 


in the Shannon information entropy and the number of principal components. Similarly, 
the manifestation of classical chaos in the statistics of quantum energy levels and the con- 
irmation of random matrix fluctuations in molecular spectra has been established in Refs. 


22 


24]. 


In the present work, we consider the quantum many-body dynamics of ten bosons con¬ 
fined in a one-dimensional harmonic trap that interact with a contact interaction potential. 
The one-dimensional regime is achieved experimentally in optical and magnetic traps with 
tight transverse confinement and frozen radial degrees of freedom. Quantum many-body 
effects are more important in reduced dimensional and interacting systems, because the 
competition between statistical properties and quantum fluctuations is enhanced in them. 
Experimentally, one-dimensional harmonically trapped quantum degenerate systems have 
been realized; see, for instance, Ref. 25|. It is hence of fundamental both theoretical and 
experimental interest to understand the time-evolution of entropy and onset of statistical 
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relaxation in the quantum many-body dynamics of the ultracold Bose gas in one dimension. 

We consider quench dynamics, where our 10-boson system is perturbed by abruptly 
switching on the two-body part of the Hamiltonian. Increasing the positive prefactor, Ao, 
of this two-body part of the Hamiltonian moves the system further and further away from 
integrability. We refer to Ao as strength of the inter-particle interactions. The process of 
statistical relaxation is then studied by analyzing the evolution of the Shannon entropy, 
occupation entropy, and the number of principal components in the time-dependent many- 
body basis set encompassing the multiconhgurational time-dependent Hartree method for 


bosons (MCTDHB) 26l,l27j (see below). Statistical relaxation is characterized by an increase 


in entropy and the vanishing of its fluctuations. MCTDHB can in principle [26, [27] and 


practice 


28 


30] provide exact solutions of the the time-dependent many-body Schrodinger 


equation (TDSE). Herein, we use the MCTDHB method for solving the TDSE at a high level 
of accuracy. Let us emphasize that exploring the solutions of TDSE for systems of interacting 
particles in an external trap potential is a fundamental problem in many-body physics. A 
key aspect of the present study with the MCTDHB approach is to enable faithful tests of 
the properties of thermalization processes in quantum many-body systems by employing 
an exact and realistic set of many-body wavefunctions as a basis for the computation of 
entropies and other relevant quantities. Let us add that the present work is restricted to the 
many-body dynamics of a pure state, i.e. the temperature of the system is absolute zero. 

We define the many-body information entropy S m *°(t) and the number of principal com¬ 
ponents Npc(t ) in terms of the time-dependent expansion coefficients of the state in the 
time-dependent MCTDHB basis. These measures are then compared with the occupation 

ues of the reduced one-body density matrix 


31 


As a key result we demonstrate that 


entropy S occu (t ) which is defined by the eigenva 
or occupation numbers, see Sec. Ill B H and Refs. 

S in f°(t) and S occu (t ) behave in a similar way: for increasing interparticle interactions, the 
many-body entropies saturate with time and their fluctuations become negligible. Hence, 


there is statistical relaxation 


20 


21| . despite the tendency of the MCTDHB description to 


adapt the basis to minimize the number of contributing coefficients. Importantly, we also 


demonstrate that S info 


t ) and S occu (t) approach the value predicted by the GOE random 


matrix ensemble 33H35I] for larger values of the interaction strength. 


To complement the results for the entropies S in ^°(t) and S occu (t) we study the time- 

,0 nge visibility 


evolution of the correlation function 


36 
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in interference experiments. For a Shannon entropy defined with a noninteracting basis set it 


was shown in Ref. 


that the interference fringes in ballistic expansion become less visible 


in the case of irregular dynamics, i.e., when the entropies become large. We investigate the 
relation of the correlation function to the production of many-body entropy quantitatively: 
we compare the time-evolution of the spatial first-order correlation function or coherence g d) 
and the many-body information entropy S in ^°{t). Our results demonstrate the close relation 


of the loss of coherence anc 
functions can be measured 


increase in many-body entropy. Since the spatial correlation 


39l. |40|], this relation can be scrutinized in experiments to test 


statistical relaxation directly. Let us recall here that testing the ETH means to verify 
that statistical relaxation is present in the system’s many-body observables fl fl Bfl. 
Measuring the expectation values of general many-body operators is however a difficult if 
not impossible task. Our results demonstrate that this problem can be circumvented, since 
statistical relaxation can be inferred from measuring the one-body correlation function. 

This paper is organized as follows. In Sec. II we give the Hamiltonian used in the present 
work and a brief introduction to the numerical many-body method MCTDHB as well as the 
quantities of interest. Sec. Ill discusses the time-evolution of many-body entropy measures 
as a function of the interaction strength to which the system is quenched. Sec. IV considers 
correlation functions and coherence in the dynamics and demonstrates a link of statistical 
relaxation to the loss of coherence. Sec. V gives conclusions of our work. 


II. METHODOLOGY 

A. Numerical method and Hamiltonian 

The evolution of N interacting bosons is governed by the TDSE, 

H^ = i—. ( 1 ) 

dt V ; 

The total Hamiltonian we consider is 

N N 

H(x!,x 2 , —,x n ) = y ^h(xj) + 9(f) y: W(xi-Xj). (2) 

i— 1 i<j = l 

Here, h(x) = T(x) + V{x) is the one-body Hamiltonian containing the external trapping 
potential V and the kinetic energy T = W(xi — Xj) is the two-body interaction of 
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particles at positions Xi,Xj, and 0(f) is the Heaviside step function of time t. The Hamilto¬ 
nian H is in dimensionless units. It is obtained by dividing the dimensionful Hamiltonian by 
— -j ^2 (m is the mass of the bosons, L is an arbitrary length scale). Since the time-evolution 
starts at t — 0, the 0(f) term in Eq. (J2]) above implements an interaction quench: the 
interactions are abruptly turned on at t = 0. We initialize the system in the ground-state 
|T(f = 0)) of the noninteracting Hamiltonian (Eq. (J2]) for t < 0). In the MCTDHB method 
which we use to solve the TDSE, Eq. (JT]) , with the Hamiltonian, Eq. (|2]), the ansatz for the 
many-body wave function is taken as a linear combination of time-dependent permanents 
with time-dependent weights, 


n 


M 

n;t) = n 

i=l 



(3) 


Here, the summation runs over all possible N con f = ( + JV _1 ) configurations, {| n;t) = 
1 77 , 1 , •••, tim] t); JT rii = N}, b\(t ) creates a boson in the ith single-particle state 0j(x, t), 
and | vac) denotes the vacuum. It is important to emphasize that in the ansatz [Eq. (|3l) ] 
both the expansion coefficients Yhi n i = N} and the orbitals {(^(x, that build 

up the permanents |n; t) are time-dependent, fully variationally optimized quantities. MCT¬ 


DHB has been established as the currently most efficient way to solve the time-depenc 


many-body problem of interacting bosons accurately and for a wide set of problems 28-30]. 


ent 


In MCTDHB(M), the vectors n = ( 771 , represent the occupations of the orbitals in 

a single configuration and preserve the total number of particles, n\ + ... + um = N. M is 
the number of single-particle functions that make up the permanents |n;t). The efficiency 
of the method comes from the variationally optimized and time-adaptive basis that makes 
the sampled Hilbert space dynamically follow the motion of the many-body dynamics. 

In the limit of M —> oo the set of permanents {[n;f)} spans the complete fV-boson 
Hilbert space and the expansion in Eq. ([3]) is exact. In practice, we have to limit the size 
of the Hilbert space in our computations. Because the permanents are time-dependent, a 
given degree of accuracy is reached with a much shorter expansion, as compared to a time- 
independent basis. This leads to a significant computational advantage over, for instance, 
exact diagonalization techniques (see explanation below). To solve the TDSE, Eq. (JTJ) , for 
the wave function l'b(f)) one needs to determine the evolution of the coefficients {Cft{t)} 
and orbitals {4>i(x, t)}^ 1 in time. Their equations of motion are derived by requiring the 
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stationarity of the action functional with respect to variations of the time-dependent coeffi¬ 
cients and the set of time-dependent orbitals. The obtained equations form a coupled set of 
nonlinear integro-differential equations 26] that we solve simultaneously with the recursive 
MCTDHB (R-MCTDHB) package 44], For reference, we also give the equations of motion 
in Appendix [XI In order to calculate eigenstates of the Hamiltonian H [Eq. (12]) ]. one uses the 
so-called improved relaxation method. By propagating the MCTDHB equations of motion 
[Eqs. (1A2|) and (IA3D ] in imaginary time for a given initial guess, excitations are exponen¬ 


tially damped and the system relaxes to the ground state. It should be noted that in the 
widespread time-dependent Gross-Pitaevskii (TDGP) theory, the many-body wave function 
is given by a single permanent |rq = N\t ); all particles reside in the same single-particle 
state (orbital) and there is consequently only a single coefficient. MCTDHB(M) contains 
the TDGP theory as the M = 1 special case. 

We stress here that MCTDHB is much more accurate than exact diagonalization meth¬ 
ods at the same dimensionality of the considered space. In exact diagonalization a time- 
independent basis is employed. In most cases, it is built from the eigenstates of a one-body 
problem. These states are not further optimized to take into account the dynamics and 
correlations in the considered system which necessarily arise due to the presence of inter¬ 
particle interactions. Thus the space and basis considered in exact diagonalization is fixed 
and not optimized, especially for the treatment of dynamics. MCTDHB on the other hand 
uses a time-adaptive many-body basis set (see Eq. (J3]), Ref. 26j, and also Appendix lAl). 


Its evolution follows from the time-dependent variational principle 


45] and is such that the 


error resulting from the truncation of many-body Hilbert space is minimized by the basis at 
any given point in time. This advantage of time-adaptivity helped to solve numerically ex¬ 


actly the time-dependent many-bod y pr ob 


em even for long-range time-dependent one-body 


and two-body potentials, see Refs. 28, [29]. These references show in a direct comparison 


of MCTDHB and exact diagonalization that the accuracy of MCTDHB for many-body dy¬ 
namics of ultracold bosons can in many cases not be reached at all by exact diagonalization 
methods, even when a very large configuration space is used. 

The Hamiltonian of bosons in one dimension is given by the Lieb-Liniger model when 
the two-body interaction is assumed to be mediated by a contact potential [46]. In the 
limit of ^ —y 0 (here n is the particle density, g is the inter-atomic coupling strength) 
fermionization occurs. This so-called Tonks-Girardeau regime is characterized by a single 
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but strongly correlated configuration [47], i.e., the expansion Eq. ([3]) would include only a 
single term. On the other extreme, the | —>■ oo limit can be captured by the TDGP mean- 
field approximation for weakly interacting bosons, i.e., a single uncorrelated configuration. 
One could hence speculate that an entropy measure related to configuration space would 
decrease as the system enters the Tonks-Girardeau regime since the eigenstates are single 
configurations. However, in our study we did not encounter such a behavior. This means that 
the dynamics in our study is far from equilibrium and adiabaticity and the Tonks-Girardeau 
states do not play much of a role. Moreover, we chose interaction strengths, for which ^ —>■ oo 
is not fulfilled sufficiently well; our investigations are therefore in the crossover between the 
two regions ^ —>■ oo and ^ —> 0 where any mean-field approach (Gross-Pitaevskii 48] or 
multi-orbital 49]) breaks down, because many different configurations are contributing. In 
this regime, as we will see in the following, many-body entropies and correlations become 
important and their features are present in the quantum dynamics. 


B. Quantities of interest 

1. Many-body entropies 

To study statistical relaxation and thermalization, we employ the measures information 
entropy S m ^°{t) and number of principal components N pc (t) defined as follows 

Nrc{t) = E«|C»(i)l 4 ’ (4) 

5 in/ ”(i) = — ElCsWtlulCWWI 2 - (5) 

n 

S info {t) and N pc {t) measure the effective number of basis states that contribute to a given 
many-body state at time t. The many-body nature of these quantities for the MCTDHB 
basis set [Eq. (J3])] can be made explicit by writing a coefficient \Cfi(t)\ 2 as an expectation 
value, 

icswr = r <t|(Si(or • ■ • ■ • • (6i,w)“"i*> • 

IL=1 ^ 

Obviously, all M creation/annihilation operators contribute to the value of every coefficient. 
For this reason, the coefficients and their distribution can also be used to directly assess the 
content of many-body entropies in the the system qualitatively. 







In the multi-orbital 


49| and Gross-Pitaevskii 48] mean-field approaches only a single 


configuration and coefficient is included and consequently S in ^°[t) = 0 and N pc (t ) = 1. 
Thus, S in f°(t) and N pc {t ) are entropies that cannot be produced in mean-field theories. The 
information entropy and number of principal components of the MCTDHB basis therefore 
are a quantitative measure for how well or not a given many-body state is captured by 
mean-field theories. Large N pc (t ) or S in f°(t) means a state contains many configurations 
and cannot be captured by mean-field methods. Small N pc (t ) or S in ^°{t) means that the 
state contains few configurations and is close to a mean-field state. 

Besides the two measures S in ^°(t) and N pc (t), we consider also the occupation entropy 


31 


, 32J defined by 


S““W = -5><W[lnnV) 


( 6 ) 


S occu (t ) is an entropy obtained from the natural occupations, i.e., the eigenvalues of the 
reduced one-body density matrix fq(f) = (see Eq. (jHJ) below and Ref. 36] for details). 
For the TDGP mean-field one has S occu (t ) = 0 always, since there is only one natural 
occupation tR = rii/N = 1 in this case. For multi-orbital mean-field theories, several occu¬ 
pation numbers can be different from 0; however, these occupations are time-independent, 
i.e., dtfii = 0. Hence, for multi-orbital mean-field theories the occupation entropy S occu (t ) 
remains constant, i.e., d t S occu (t) = 0. In the time-evolution of a many-body state, the prox¬ 
imity of the value of S occu (t ) to 0 is a measure of how well the state can be described by 
the TDGP mean-field. The magnitude of the fluctuations in S occu (t) indicates how well the 
state could be described by a multi-orbital mean-field approach. Since there are non-mean- 
field states with many configurations that also fulfill d t S occu (t) = 0, it has to be checked if 
additionally S in ^°{t) rs 0 holds in order to conclude that a state with constant occupation 
entropy is indeed of single-configuration (i.e., mean-field) type. 

For a GOE of random matrices, N^° E = D/3 and Sqq° e = In 0.48D holds, where D x D 
is the dimension of the random matrices. To obtain the GOE estimate for the information 
entropy S info (t ), we set D equal to the number N con f of time-dependent many-body states 
in MCTDHB. For the occupation entropy, the GOE analog is obtained by setting n* = jj 
for all i in Eq. ([6]) . Consequently, in the case of S occu (t), D is set equal to the number of 
orbitals M in the MCTDHB treatment. We obtain Sqq u e = — {j!) (if) = l n (^ /j 0 i n 

this case. 
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2. Correlation functions and coherence 


In the present work, we investigate the normalized first-order correlation function 
defined as 36, 37] 


(i)/ / p (1) (xi|xi;i) 

g y ’{x-tjXi]t) = , = 


(7) 


where p is the diagonal part of the one-body density matrix pd) gi ve n by 


p (1) Oi|xi;t) = 

N J if*(x [,x 2 ,x N ; t)ip(x 1 ,x 2 ,..., x N ; t)dx 2 dx 3 ... dx N . (8) 

The normalized spatial first order correlation function pd) quantihes the degree of first 
order coherence. Note that Ig^^x^, xp, £)| < 1 means that the visibility of interference 
fringes in interference experiments will be less than 100%; this case is referred to as loss of 
coherence. \g ( ' 1 \x' 1 , xp,t)\ = 1 corresponds to maximal fringe visibility of the interference 
pattern in interference experiments and is referred to as full coherence. The strength of 
the interaction between the particles affects the correlations: the stronger the inter-particle 
repulsion, the stronger the loss of coherence. 


III. THE PRODUCTION OF MANY-BODY ENTROPY IN INTERACTION QUENCH 
DYNAMICS 

Our present calculations are performed for one spatial dimension and consider N = 10 
repulsively interacting bosons in an external harmonic trap. We use a contact interaction 
W(xi — Xj) = Ao 5(xi — Xj ) and the external trap V(xi) = ^xj in Eq. (J2]). The dimensionless 
strength Ao of the inter-particle repulsion is varied in our investigation. We restrict the 
number of orbitals to M = 6, yielding a total of D = N con f = ( iV+ ^ -1 ) = 3003 permanents. 

A key aspect of many-body quantum chaos is how the many-body wave function spreads 
across the available basis functions with an increase in the inter-particle interaction strength. 

As mentioned, this spreading is counteracted by the time-adaptivity of the basis set in 
MCTDHB. 

For an assessment of the presence of statistical relaxation, it is instructive to first visualize 
the time-evolution of the coefficients directly. In Figs. 1 and 2 we plot \Ca(t)\ 2 as a 
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function of the index n of the basis states | n; t). See Ref. 50] for the details on how to obtain 
the index n from the vectors n = {n i, Fig. 1 is for the interaction strength Ao = 0.5 

and for times t = 0.2, 0.5,1.0, and 2.0. Note that the permanents | n\t) are not eigenstates 
of the Hamiltonian for t > 0. However, every eigenstate of the interacting problem can be 
represented as a pattern of coefficients that contribute significantly. Initially, at t — 0, only 
a single coefficient is nonzero, because the system was relaxed to the noninteracting ground 
state. As depicted in Figs. la)-d), the number of significantly contributing coefficients grows 
with time due to the inter-particle interactions, but remains only a small portion of the 
available basis states for M — 6. This is a consequence of the interactions being relatively 
weak. If the number of non-zero elements of {Cfi(t)\ JAn* = N} is only a small portion 
of N con f = 3003, we refer to the respective state as localized. Localized states are rather 
close to a mean-field description for which only a single coefficient would contribute. The 
results in Fig. 2 are for a larger value of the interaction strength (Ao = 10.0). It is seen from 
Fig. 2 that there is a substantially larger amount of non-zero coefficients as compared to 
the smaller interaction strength Ao = 0.5 (Fig. 1). Since the contributing coefficients spread 
over almost the whole available space, we refer to such a state as delocalized. Delocalized 
states cannot be captured by mean-held descriptions. From Figs. 1 and 2 we conclude that 
increasing the interaction strength makes the states emerging in the time-evolution change 
character from localized to delocalized. 

To further quantify the average number of time-dependent basis states that make up the 
wave function, we plot exp(S' m ^°(t)) and N pc (t ) for the small interaction strength Ao = 0.5 
in Fig. 3, and for the large interaction strength Ao = 10.0 in Fig. 4. For Ao = 0.5, the 
inter-atomic correlations are small and large fluctuations in both N pc (t ) and exp(5' m ^°(t)) 
are seen (Fig. 3). These fluctuations are around comparatively small averages (compare 
Figs. 3 and 4). We conclude that for localized states, there are only a few contributing 
components. For larger interaction strength (Fig. 4), the fluctuations in exp (S in *°(t)) and 
N pc (t) decrease. This is analog to the GOE behavior and in line with computations employ¬ 


ing a time-independent basis set 


M 


35 


511 ]. An important observation from Fig. 4 is that 


the entropies saturate for large interaction strength. We can conclude from Figs. 3 and 4 
that increasing Ao makes the fluctuations of the number of principal components N pc (t ) and 
exp(5' m ^°(t)) decrease and results in the emergence of a saturation of these quantities. 

To shed more light on this behavior, we continue by comparing the many-body entropies 
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S info {t ) and S occu (t). 

Fig. 5 shows S m f°(t) and S occu (t ) (Eq. (jSJ) and (EJ), respectively) for the interaction 
strengths A 0 = 0.5,10.0, and 15.0. As a first observation, we find that S m f°(t ) and S occu (t ) 
show a similar overall behavior. We discuss first the information entropy S in ^°(t) in detail 
and subsequently S occu (t). 

For smaller interactions, Ao = 0.5, the increase in S in ^°{t) is almost linear and reaches 
the value S in ^° ~ 3.2 at time t — 6 [Fig. 5a)]. The GOE value S™q E for our present 
calculation with 10 bosons, 6 orbitals, and D = N con f = 3003 is Sqq° e = ln(0A8D) = 7.273. 
Since we observe S ini °(t) < Sqq° e , we infer that the system remains relatively ordered. 
Nevertheless, the value of is far enough from 0 to conclude that the state can no 

longer appropriately be described by mean-field methods. For larger interaction strengths, 
S info (t) shows quick saturation [Fig. 5b) and c)]. For A 0 = 10.0, saturation emerges close to 
gmfo _ g 4 g anc [ f or _ 15 g gq gmfo _ y.17. The values of S m f°(t) are approaching the 
GOE value Sqq° e = 7.23. The discrepancy is because the interaction strength, albeit being 
comparatively large, is still finite. As far as time-independent basis sets are concerned this 
may be attributed to the operation of the so-called embedded GOE of random matrices in 
interacting particle systems (Ref. [52]). The theory of embedded GOEs can be applied to 
systems with lower-body rank operators such as the one- plus two-body Hamiltonian used 
in the present study. It predicts that entropy measures will be close to but not identical 


to the GOE value |33[]. Importantly, the saturation of the many-body information entropy 
S in f°(t) close to the GOE value Sqq° e can be seen as a hallmark of statistical relaxation. 
This saturation demonstrates that statistical relaxation gradually overcomes the tendency 
of the time-dependent variational principle to minimize the spread of the coefficients by 
optimization of the MCTDHB basis set as the inter-particle repulsion Ao increases. Let us 
now consider the occupation entropy S occu (t ) in more detail. 

In Fig. 5, the results for S occu (t ) are shown for the interaction strengths Ao = 0.5,10.0, 
and 15.0. At t — 0, all bosons are in the lowest orbital: n\ = N and rii = 0 for i = 2, 3...M. 
Then ni = 1 and all other h* = 0 [cf. Eq. (jBJ)]. Therefore S occu (t = 0) = — ln(l) = 0 as in 
the case of the information entropy S in ^°[t = 0). The system is, hence, in a mean-field state 
at t — 0 (cf. Sec. Ill B If) . At later times, however, the bosons start to distribute themselves 
in all M orbitals. For all interactions [Figs. 5a),b), and c)] S occu (t ) is saturating after 
an initial adjustment to the quench. The state is, however, not of (multi-orbital) mean- 
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field character, since S in ^°(t) > 0 indicates that several configurations are contributing 
to the dynamics. For small interactions [Fig. 5a)], S occu {t) saturates far below the value 
of the GOE. The state remains rather localized and the dynamics can be considered as 
regular and we refer to the state as ordered. For larger values of the interaction strength 
(A 0 = 10.0 and A 0 = 15.0 in Figs. 5b) and 5c), respectively), when the wave function becomes 
fully delocalized with time, all bosons are on average roughly equally distributed in the M 
orbitals. The situation resembles the GOE and S occu (t ) saturates close to its GOE estimate 
Sqqe■ I 11 our calculations with M = 6 and N = 10, Sqq u e = In(M) = ln(6) = 1.79, as 
discussed in Sec. Ill BT1 This convergence of S occu {t) to Sqq u e for increasing Ao demonstrates 
the presence of statistical relaxation and is similar to the many-body information entropy 
S mfo (t) discussed in the previous paragraph. As also argued for S m ^°(t), the value at which 
S occu (t ) saturates is smaller than the GOE value Sqq u e due to the applicability of the theory 
on embedded GOEs for large but hnite interactions. One may speculate that the chaotic or 
irregular dynamics which we find are a many-body analog of the wave chaos found for the 
TDGP case in Ref. 32]. 

Fig. 6 shows the information entropy S in ^°[t) for much longer times than Fig. 5. For 
interaction strength Ao = 0.5, S in ^°{t) does not reach saturation and exhibits strong fluctu¬ 
ations for a long time [Fig. 6a)]. The system may eventually reach some equilibrium state 


which is different from a thermal state and may, according to Ref. 


531] . be described in terms 


of a generalized Gibbs ensemble. For Ao = 15 it is seen from Fig. 6b) that there is a plateau 
at small times (t ~ 0.5, see Fig. 6b) inset). The relaxation to the GOE value Sqq E happens 
on a much longer time scale (t ~ 10) and is non-transient. We speculate that the first stage 
of relaxation on a much shorter time scale may be a signature of so-called prethermalization. 
Prethermalization is characterized by a rapid relaxation of only some observables (not all, 
as in the case of the ETH) to their equilibrium values [54] and has recently been observed 
experimentally for a degenerate one dimensional Bose gas 55]. 


IV. RELATION OF CORRELATION FUNCTION AND COHERENCE TO MANY- 
BODY ENTROPY PRODUCTION 

In this section we discuss the time-evolution of the first-order correlation function cpb 
(see Sec. Ill B 21 for a definition) for our system. The correlation function quantifies the 
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coherence which can be experimentally determined in interference experiments. We would 
like to put forward a strong link between the dynamics in the coherence and in the many- 
body measures of entropy defined in Sec. Ill B II and analyzed in the previous Section. This 
connection can be used to quantify the many-body entropy in the system by measuring the 
correlation function. 

Fig. 7 presents aq; t) | 2 as a function of its two spatial variables for two interaction 

strengths, Ao = 0.5 and Ao = 10 at various times t. For weak interactions, i.e., Ao = 0.5, 
\g ( ' 1 \x' 1 , x^t) : remains close to unity for all for a comparatively long time. The 

system remains coherent throughout its time-evolution. This observation is in-sync with 
the small spreading of the states coefficient distribution (see Fig. 1) for small interactions 
A 0 = 0.5. Turning to stronger interaction strength (A 0 = 10.0), after a sufficiently long 
time (t ~ 10) the correlation function is unity almost only along the diagonal (x\ = xf). 
Away from the diagonal (x\ ^ aq) the correlation function \g^ l \x'^ aq; t)| 2 is close to 0. 
Hence, for the stronger interactions Ao = 10.0, the coherence of the system is lost with time. 
That strong inter-particle repulsion leads to an almost complete loss of coherence agrees 
with our previous observation of statistical relaxation, saturation of entropies S m ^°{t) —> 
Sgoe an d S occu (t ) —> Sqq u e , and delocalization of the MCTDHB coefficients distribution 
(compare Figs. 1,5, and 7). We conclude that the production of many-body information 
entropy S in ^°(t) and occupation entropy S occu (t ) entails an intensified loss of coherence. 
This connection can be exploited to measure the many-body entropy of a system, because 
the one-body correlation function can be determined through interference experiments. 
Let us stress again here that the larger the disorder [measured by the entropies S in f°(t) and 
S occu (t)\ in a quantum mechanical system the less it can be described by a product of a 
single complex valued function. Since the ansatz of the TDGP uses a single complex-valued 
function, our finding implies that the TDGP can generally not adequately describe processes 
such as statistical relaxation or thermalization which are characterized by a loss of coherence 
and an in-sync increase and saturation of the entropies S m ^° and S occu . 


V. CONCLUSIONS 

We have studied many-body entropy production, statistical relaxation and coherence of 
parabolically trapped interacting bosons for an interaction quench by highly accurate MCT- 
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DHB computations going beyond the scope of commonly applied mean-field approaches. The 
full time-dependent solution of the many-body problem with MCTDHB allowed us to de¬ 
fine and compute new many-body entropy measures. We analyzed the information entropy, 
number of principal components, and occupation entropy of the time-dependent MCTDHB 
basis set as a function of time. We have shown that an increase in the many-body entropy 
measures is linked to a loss of coherence in the dynamics (see Figs. 5-7). For larger values 
of the interaction strength, we find irregular dynamics and statistical relaxation despite the 
tendency of the time-adaptive many-body basis set of MCTDHB to minimize the number 
of contributing expansion coefficients and therewith the entropies related to them. All the 
entropy measures are in mutual agreement (see Figs. 3-6). 

For larger inter-particle repulsion, the expansion coefficients \Cft(t)\ 2 delocalize more 
strongly (Fig. 2) and we observe a quick initial saturation of many-body entropies (see 
Figs. 5 and 6) possibly related to prethermalization. A saturation of the many-body en¬ 
tropy to the GOE value follows on a longer time scale. It is important to stress here that this 
saturation is not as trivial as one might be tempted to assume: the basis set in MCTDHB 
is explicitly time-dependent and it is optimized such that it minimizes the portion of signif¬ 
icant coefficients in the expansion. Nevertheless, the many-body entropy of the coefficients 
and the eigenvalues of the reduced one-body density matrix approach the GOE values and 
statistical relaxation prevails. 

By studying the time-dependence of the Erst order correlation we demonstrate a strong 
link between the dynamics of entropy and the dynamics of coherence. Our present work ex¬ 
emplifies that large production in many-body entropy causes an intensified loss of coherence 
(Figs. 5-7). This loss of coherence constitutes an independent signature of the existence of 
statistical relaxation, allowing to study it from another perspective and, most importantly, 
to measure it in experiments. 


Further investigations are needed to test the ETH and to assess the generality of statistical 
relaxation as well as its relation to chaos, also for larger portions of many-body Hilbert space. 
In this respect it is interesting to understand if the found emergence of statistical relaxation 
is a many-body analog of the wave chaos in the TDGP equation reported in Ref. [32]. Finally, 
a straightforward continuation of this work would be to hnd its connection to investigate 


the possible connection the recently observed prethermalization phenomenon 


55]. 
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Appendix A: The MCTDHB equations of motion 


In the following, the equations of motion of MCTDHB are given and their derivation is 


sketched. For details, see Ref. 


26]. The action of the TDSE, 


S= dt l(^\H-id t \ty)+ ^2^ ij (t)((<f) i \<f) j ) - Sij)] , 


(Al) 


V 


is demanded to be stationary with respect to the variation of the time-dependent orbitals 
as well as with respect to the variation of the time-dependent coefficients 
{Cn(t)\Yhi n i = -A}- The orthonormality of the orbitals {</>*,(£; f)}^ is enforced by the 
Lagrange multipliers /i^ (t) in S. From the stationarity of the action S the equations of 
motion of the orbitals, 


M 


idt^j^X.t) = P ( h4j{x,t) + Y {pityjkPksqlWWslfatWqfat) ), 

\ k,s,q,l =1 

M 

P = 1 - Y 


(A2) 


j '=i 


as well as of the coefficients, 


'H(t)C(t) = id t C(t)-, = (rh]t\H - id t \m';t), (A3) 
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are derived. These equations flA2j) . flA3j) form the core of MCTDHB. Without loss of gen¬ 
erality, equations f1A2[) , flA3[) are given here for one spatial dimension under the constraint 
(dt4>k\4>j) — 0 Vk,j, compare Ref. 26]. The following notations were invoked for the respec¬ 
tive matrix elements and operators: 

Pk q = mblb q \*)> 

Pksqi = (ty\blblb q bi\y), 

W s i(x-t)= I dx'(j>*(x,t)W(x,x , -,t)<i>i(x , ,t), 

h =-\b +v(x) - 

Finally, the shorthand notation C(t) collects all the time-dependent coefficients {C*j(i); n\ Vb n* 
N} in a vector, employing an enumeration scheme documented in Ref. [50]]. The MCTDHB 
equations of motion form a coupled set of nonlinear integro-differential equations because 
the evaluation of the matrix elements Pk q , Pksqi in the orbitals equations (1A2|) depends on 
the coefficients C(t), and the evaluation of the coefficients equation (1 A - > I) depends on the 


matrix elements of the Hamiltonian with the current set of orbitals, hk q = (4>k\h\<j> q ) an d 
Wksqi = J dx f dx / </)k(x,t)(/)*(x,t)W(x,x';t)(/>i(x , ,t)(/>g(x,t). The equations (1A2jl . (1A3|1 can 
be solved efficiently and self-consistently with the R-MCTDHB pa ckag e 44], For further 
details on the derivation and properties of these equations, see Ref. [261. 
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FIG. 1. Time-evolution of the distribution of the magnitude of the coefficients {|C^(t)| 2 } for 
a weak interaction, Ao = 0.5. Panels a)-d) show the magnitude of the coefficients at times 
t = 0.2,0.5,1.0, 2.0, respectively. The index n is computed from the vector ft using the mapping 
described in Ref. 50j. With increasing time, more coefficients in the expansion become significant, 
but the spread is far from the whole available space spanned by the D = N con f = 3003 configu¬ 
rations. The state stays rather localized throughout the quench dynamics. The coefficients with 
n > 70 are smaller than 10 -4 and not plotted therefore. All quantities are dimensionless. 
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FIG. 2. Time-evolution of the distribution of the magnitude of the coefficients {|C^(t)| 2 } for 
a strong interaction, Ao = 10.0. Panels a)-d) show the magnitude of the coefficients at times 
t = 0.2,0.5,1.0, 2.0, respectively. The index n is computed from the vector ft using the mapping 
described in Ref. 50|. With increasing time, more coefficients in the expansion become significant; 
the coefficients explore almost the whole available space spanned by the D = N con f = 3003 
configurations. The state rapidly becomes rather delocalized throughout the quench dynamics as 
compared to the localized case in Fig. 1. All quantities are dimensionless. 
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FIG. 3. Time-evolution of the number of principal components N pc (t) and the exponential of 
information entropy exp (S m ^°(t)) for the weak interaction Ao = 0.5. Both exp(5' m ^°(t)) (green, 
upper line) and N pc (t ) (red, lower line) exhibit the same overall behavior. In analogy to the GOE 
for small interactions, large fluctuations emerge in both quantities due to the absence of strong 
correlations between the particles. See text for further discussion. All quantities are dimensionless. 
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FIG. 4. Time-evolution of the number of principal components N pc (t) and the exponential of 
information entropy exp (S in f°(t)) for stronger interactions, Ao = 10.0. As in Fig. 3 for weaker 
interactions, both exp(S in f°(t)) (green, upper line) and N pc (t ) (red, lower line) exhibit the same 
overall behavior. The thin black dashed lines are provided to guide the eye and to estimate the 
magnitude of fluctuations. In analogy to the GOE for larger interactions, the relative fluctua¬ 
tions are quenched due to the presence of strong correlations as compared to the case of smaller 
interactions values in Fig. 3. See text for further discussion. All quantities are dimensionless. 
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FIG. 5. Dynamics of many-body Shannon entropy S in f°(t ) and occupation entropy S occu (t) for 
different inter-particle interaction strengths. Statistical relaxation causes the convergence of both 
quantities, S m f°(t ) (red, upper lines) and S occu (t) (green, lower lines), to the values S 1 qq° e and Sqq % e 
(horizontal black dashed lines), respectively, as the interaction strength increases from panels a) 
to c). Implications are discussed in the text, curves are smoothened for clarity of presentation, all 
quantities shown are dimensionless. 
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FIG. 6. Long-time-evolution of many-body Shannon entropy S in f°(t) for weak Ao = 0.5 and 
stronger Ao = 15.0 interactions. To verify that the saturation of the many-body entropy measures 
is indeed non-transient, their time-dependence for longer times is depicted. Panel a) shows S in ^°(t) 
for Ao = 0.5 and panel b) for Ao = 15.0. The behavior for smaller times is presented in the insets. 
The saturation happens on a much faster time-scale when the interaction strength is large with 
Ao = 15.0 as compared to the situation when the interaction strength is small with Ao = 0.5. 
See text for details, curves are smoothened for clarity of presentation, all quantities shown are 
dimensionless. 


23 


















FIG. 7. Coherence in the quench dynamics measured with the correlation function \g( 1 \x' 1 , x±', t)\ 2 
for weak Ao = 0.5 and strong Ao = 10.0 interaction strengths. The left column depicts [g^| 2 for 
t = 0.1,5.0,10.0 with Ao = 0.5, respectively. The right column shows [g^) 2 for the same times 
but for Ao = 10.0. The states with localized coefficient distributions and small entropies [compare 
Fig. 1 and Fig. 5a)] are also closer to being coherent ( ~ 1 ) than the states with large spread 

in coefficient distribution and large entropies [compare Fig. 2 and Figs. 5b) and 5c)]. In the case of 
spread-out coefficients and large many-body entropies [t/ 1 ^ 2 ~ 0 holds almost everywhere but for 
the diagonal \g ( ' 1 \x,x]t)\ 2 . Entropy production and loss of coherence hence go in-sync. See text 
for further discussion, all quantities shown are dimensionless. 
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